/*
 * EquilibriumSolver.h
 *
 *  Created on: 26 Jul 2011
 *      Author: Allan
 */

#ifndef EQUILIBRIUMSOLVER_H_
#define EQUILIBRIUMSOLVER_H_

// C++ includes
#include <vector>
#include <functional>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "EquilibriumCondition.h"
#include "EquilibriumProblem.h"
#include "../Multiphase.h"
#include "../MultirootSolver.h"
#include "../CanonicalReactionSystem.h"

class EquilibriumSolver 
{
public:
	struct Diagnostic
	{ unsigned iterations; bool converged; };
	
public:
	EquilibriumSolver();
	
	EquilibriumSolver(const Multiphase& multiphase, const CanonicalReactionSystem& reactions);
	
	void
	Initialiaze(const Multiphase& multiphase, const CanonicalReactionSystem& reactions);
	
	void
	SetEquilibriumConditions(const vector<EquilibriumCondition>& eConditions);
	
	void
	SetSolverOptions(double tolerance, unsigned maxiterations);
	
	const Diagnostic
	Solve(double T, double P, VectorXd& n);
	
	const EquilibriumCondition
	FixedMass(const string& species, double mass);
	
	const EquilibriumCondition
	FixedTotalMass(const string& species, double totalmass);
	
	const EquilibriumCondition
	FixedActivity(const string& species, double activity);
	
	const EquilibriumCondition
	FixedAcidity(double pH);
	
	const EquilibriumCondition
	FixedPartialPressure(const string& gas, double Pg);
	
	const EquilibriumCondition
	ChargeBalanced();
	
private:
	/// The multiphase object
	Multiphase multiphase;
	
	/// The system of equilibrium-controlled reactions in the canonical representation
	CanonicalReactionSystem reactions;
	
	/// The equilibrium problem object which overloads the methods Function and Jacobian from MultirootProblem 
	EquilibriumProblem problem;
	
	/// The multiroot solver object used to solve the non-linear system of equations
	MultirootSolver solver;
	
	/// The number of primary (Nj) and equilibrium (Ne) species
	unsigned Nj, Ne;
	
	/// The indexes of the primary (jIndexes) and equilibrium (eIndexes) species
	vector<Index> jIndexes, eIndexes;
	
	/// The primary (jSpecies) and equilibrium (eSpecies) species 
	vector<string> jSpecies, eSpecies;
	
	/// The stoichiometric matrix of the equilibrium-controlled reactions
	MatrixXd vej;
	
	/// The tolerance for the iterative procedure
	double tolerance;
	
	/// The maximum number of iterations for the iterative procedure
	unsigned maxiterations;
	
private:
	const EquilibriumCondition
	FixedMassHelper(const string& species, const FunctionTPna& h);
	
	void
	UnderRelax(VectorXd& n, const VectorXd& nOld) const;
};

#endif /* EQUILIBRIUMSOLVER_H_ */
